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Abstract 


For space missions of future, completely autonomous robotic machines will 
be required to free the astronauts from routine chores of equipment maintenance, 
servicing of faulty systmes, etc. and to extend human capabilities in hazardous 
environments full of cosmic and other harmful radiations. In places of high 
radiation and uncontrollable ambient illuminations, T.V. camera based vision 
systems cannot work effectively. However, a vision system utilizing directly 
measured range information with o ^ime of flight laser rangefinder, can 
successfully operate under these environments. Such a system will be independent 
of proper illumination conditions and the interf^'^ing effeci-s of intense 
radiation of all kinds will be eliminated by the tuned input of the laser 
instrument. 

Especially important is the capability of the laser based vision system 
to afford a 3-dimensional description of the environment. Known objects can 
be recognized and their correct spatial locations and orientation can be deter- 
mined with respect to the sensing instrument, and the robot, which can thus 
perform its tasks in a three dimensional space. Vision systems using 2 dimen- 
images as the primary input have not met much success in 3-dimensional scen^ 
description because of the difficult problems associated with indirect range 
calculation methods. 

Processing the range data according to certain decision, stochastic 
estimation and heuristic schemes, the laser based vision system will recognize 
known objects and thus provide sufficient information to the robot's control 
system which can develop strategies for various objectives. 
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I. INTRODUCTION 

Various schemes for data acquisition, edge detection, plane surface extrac 
tion by clustering, and noise removal from the raw input data were discussed 
earlier in [5]. This report extends the previous work of surface description. 

A procedure is developed for defining the surfaces by regression planes passing 
through the three dimensional points clustered by the Hierarchical Clustering 
Scheme using orthogonal surface slopes along orthogonal directions. 

Continuous inner edges are calculated by calculating the possible lines of 
intersections of the planes found from the surface description procedure. This 
procedure for calculation of edges of objects can be applied only if the two 
surfaces forming it are visible from the location of the laser sensor. The 
edges detected by the Rapid Estimation Scheme, [1], are in the form of scat- 
tered points on both sides of the edges and some of the edges especially the 
inner edges, may be missed, therefore a more elaborate procedure than the 
Rapid Estimation Scheme alone was required. Thus the inner edges are calcu- 
lated here from the intersection of best fitted planes. The calculated edges 
are verified in the case of inner edges also detected by the Rapid Estimation 
Scheme. 

For the description of a scene in terms of meaningful objects, heuristic 
rules are developed for combining planes and edges found above into consistent 
subjects so as to define meaningful objects, known or unknown. Based on the 
convexity properties of plane faceted objects different planes are linked to- 
gether if there is a convex edge between them. Similarly, half edges with one 
end on a virtual vertex are linked together if the unit vectors along them are 
col inear. Heuristic rules are also developed for linking half planes resulting 
from an occlusion. 
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Noise removal from the raw data and estimation of range slopes and, in 
turn, actual surface slopes is achieved by the applicator of a two-dimensional 
smoothing algorithm, [2]. The smoothing algorithm used three different 
smoothness measures which were reduced to a quadratic form in the state vector 
consisting of the function values, two first partial and the first mixed partial 
derivatives at each node point of a uniform two-dimensional grid. The weighting 
matrix was, however, derived by a complex reduction procedure. 

To simplify the derivation of the results, the weighting matrix form of 
the smoothing integrals was avoided and the two-dimensional smoothing problem 
was formulated as an optimal control problem by incorporating a control term in 
the recursive model of the spline functions. Essentially the same results were 
obtained for the smoothing algorithms as in, [2], thus providing mathematical 
rigor and verification. 

To employ different spline function with certain interesting smoothing or 
edge enhancement properties a general method for the calculation cf bases of 
different spline functions was developed. By using this approach, it is hoped 
that the derivation of different spline functions will be simplified. 
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II. SUMMARY OF PROGRESS 

The first part of this report presents a scene analysis/object recognition 
system that utilizes directly measured rcinge information from the scene of 
interest and provides a three-dimensional description of the scene in terms of 
known and unknown objects in their true relative positions and orientations. 

The range to a large number of points in the scene is measured by using a time 
of flight laser rangefinder in which range is obtained from the phase difference 
between the transmitted and reflected laser beams. The rangefinding system, 
in fact, measures the spherical coordinates of various points on the surfaces 
of objects by scanning in the elevation and azimuth directions in a suitably 
defined coordinate system. Figure 1. The range readings are corrupted with 
measurement noise. 

By uniform increments of the two pointng angles of the laser beam over 
suitable intervals, a matrix of range values is obtained, the rows and columns 
of which are indexed with constant values of elevation and azimuthal angles. 
Range slope in the inpath, i.e., by the variation of elevation aigle only, and 
in the crosspath direction, i.e., by the variation of azimuthal angle only, is 
defined as the first difference of ran<e divided by the corresponding increment 
in pointing angles Band 9 respectively. The range data is processed along 
columns and rows to detect the presence of edges in the scanned scene which 
manifest themselves as sudden changes in range or range slope between adjacent 
points. The edges are thus extracted by a Rapid State Estimation Scheme (RES), 
[1], that uses Kalman filters and decision trees. 

The observed range data is uniformly spaced in elevation and azimuthal 
angles, but if projected onto any of the planes xy, yz, or zx shown in Figure 1, 
the regularity of data spacing is lost. Thus the range and two range slopes 
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are smoothed in the original spherical coordinates before any scene analysis 
procedure is applied. This is accomplished by using a two-dimensional smoothing 
scheme [2], in the region free of the edges of objects detected above. 

The smoothed range and two range slopes are transformed to the actual slope of 
the surfaces of the objects at their observed data points. 

Using the fact that the slopes of a plane along any two orthogonal directions 
and with respect to a reference plane are constant at all of its points, the 
data points belonging to plane surfaces can be separated. Using the two slopes 
along orthogonal directions as the two components of a feature vector, the plane 
surfaces are extracted by using a hierarchical agglomerative clustering procedure, 
[ 3 ], [ 4 ]. 

The continuous edges of a plane faceted objects are then obtained as tie 
intersections of their plane surfaces found above. Using heuristic rules, tha 
edges and planes belonging to individual objects are grouped in an ordered and 
consistent manner so as to facilitate the object reconstruction. 

The characteristics of an object which are important for this purpose are 
found to be the convexity and colinearity of edges. Methods for determining if 
an edge is convex or concave and if two edges are col inear are developed. These 
methods are consistent, regardless of the viewers point of observation. If parts 
of the image of the workscene in the input information exhibit certain proper- 
ties, then these parts will be grouped together as an object. 

For recognition of known objects, three-dimensional models in the form of 
ordered li«ts of features of objects are stored in the computer memory. The 
essential contents of these models are the specification of number and lengths 
of edges meeting at a vertex, the spatial angles between these edges, the edges 
bounding a face, etc. This stored information can then be utilized to recognize 
known objects and determine their positions and orientations. 
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The second part of the report deals with the problems of two-dimensional 
filtering and smoothing of data by using spline functions. A number of ap- 
proaches to the 2-D data processing problem have been proposed in the literature 
recently. These algorithms typically fall into two classes, those based on 
original image model [6-8], and those found by minimizing objective function 
without explicit model [2], [9-12]. The approach taken here is closer to the 

technique of numerical spline analysis, that is, no prior image model is 

assumed. However, the focus here is on achieving a recursive algorithm which 

can be Implemented on-line. 

The purpose here is twofold: a) along the same line as developed in [2], 

the results were extended to a wider class of possible signals; b) the funda- 
mental principles and limitations on recursive smoothing splineL are explored, 
which were not discussed in the previous research. 

The state space notation of signal is used for it not only keeps the nota- 
tion simple, compact, and unique, but also provides the crucial point in the 

recursive processing. Another factor involved in the formulation of recursive 
computation is the proper selection of the objective function. 

In particular, it is hoped that the expositions will shed some light on the 
peculiar probmes of recursive algorithms for 2-D smoothing splines, which do 
not appear in 1-D data processing. 

In order to get better data fitting, p.H.p.'s of continuity greater than 
the first derivatives are desirable. However, the derivation of high order p.H.p. 
using direct method is quite complex in general. A method using cardinal bases, 
[13], is also described which results in a systematic, simple, and straight- 
forward calculation procedure. 
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III. DETAILED PROGRESS REVIEW 

The progress made during the past year In the ongoing research Is described 

here. There are essentially two major phases of research, the first one dealing 

with the three-dimensional scene analysis probltm and the other Is on the 

processing of two-dimensional signals. 

TASK A. Scene Analysis Based on Laser Rangefinder Data 
Processing for Space Robotics. 

This part of the report deals with the development of a possible vision/ 
environment sensing system that can be used by autonomous robotic systems of 
future, especially In space probes where proper illumination for T.V. cam<»ra 
based vision systems will be difficult. 


1 . Feature Extraction and Surface Description 

Application of the two-dimensional smoothing algorithm [ 2 ], using spline 

functions, in the regions between the discrete edges, detected by the Rapid 

Estimation Scheme, [ 1], yields smoothed values of range r.. and the two range 
^'"ij "'“ij 

slopes 3^ & af'd the measured pointing angles and at an observed 

point P^j, Figure 2a, the actual slope of the surface at P^j with respect to 
some coordinate planes xy, yz, zx of the laser rangefinder and the three pairs 
of reference directions employed are as shows by 


a. Segmentation of Range Data: Recursive Hierarchical Clustering 

The Idea behind the transformation from range slopes to the actual surface 
slopes of the objects Is that the data points belonging to a plane surface will 
show some slopes along two fixed orthogcr.^.l directions and thus can be separated 
from the rest of the points. Hence, for all of the observed data points, two 




slopes along any two fixed directions w»r«t a reference plane cnan serve to 
separate the total observed data into plane surfaces. However, there can arise 
a serious problem due to some slopes becoming singular or achieving arbitrarily 
large magnitudes. This will happen when the corresponding unknown plane has a 
high inclination w*r-t the reference plane (Figure 2b). 

Thus for avoiding the problem of singularities in the transformed surface 
slopes it was decided to calculate the same along six different directions, 
Figure 2a. For a particular observed point only a pair of surface slopes w*r.t 
one of the coordinate planes is calculated. Possible pairs of slopes are: 



As can be noted, only three of the above quantities are unique. Denoting the 

ar.i 3r.. 

range slopes by and by three of the above surface slopes 
i J i J 

may be written in terms of the measured and smoothed quantities: 



d£ ^ 

-^tanBlV - r 

(2) 

dy 

[cose - r tan 6 sine /VC] V - r tanS cose 

$2 A 


-[tane + r/V]V^ 
r cose + sine (1 - r tan 6/V)Vc 

(3) 


^ . 

[tane + r/V^] V + r tane tans 

(4) 

dx 

[1 - r tane /VC] v - r tan 6 


where each of the above expressions connect all the quantities at a single point. 
The remaining three surface slopes In Equation 1 may be obtained as the recip- 
rocals of these quantities. 

Thus, with each observed point of the scene, a two-dimensional slope vector 
is associated indexed by the particular plane w-r-t which the two slopes were 
calculated. Processing the whole data in this manner, 6 two-dimensional arrays 
corresponding to six different slopes are obtained. These arrays contain blank 
entries for points at which the particular slope was not calculated. 
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The observed data points are separated into different groups for dif-Perent 
planes surfaces by using a recursive hierarchical clustering scheme reported 
earlier, [5]. 

b. Surface Estimation: Least Squares Plane Fitting 

The output of the recursive clustering algorithm [5], consists of three 
different sets of clusters, depending on which pair of surface slopes was used 
to group the data points. For a cluster of three-dimensional points using 
surface slopes ^ and 0 a linear relation z = a + bx + cy, (a plane) is 
supposed to exist between their rectangular coordinates (x,y,z). Given 'n' 
number of points, the best values of the parameters a,b, and c are obtained by 

Q 2 2 

minimizing the squared error 1 (z. - z) = E (z. - a - bx. - cy.) (5) 

i=l ^ i=l ^ ^ ^ 

Differentiating this error w-r-t a,b, and c, three euqations are obtained: 


If = (z, - a - bx, - cy,) (-2) (6) 

II = (z, - a - bx, - cy,) (-2x,) (!) 

If = (z, - a - bx, - cy,) (-2y,) (8) 


Setting all derivatives equal to zero for a = S, b = B, and c = c, the so-called 
normal equations of the least squares plane are: 

Z z. = na + S Z x^. + c Z y. (9) 

Z x^.z^. = a Z X 2 + b Z x^^ + c Z x^.y^. (10) 

r = a Z y . + b Z x.y. + c Z y^^ (11 ) 

where all summations are from i=l,....,n. Solution of this set of linear 
equations defines the best fitting plane as the triple (a, 6, c). 
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c. Continuous Edges: Calculation and Verification 

The edge points detected by the Rapid Estimates Scheme, [1 ], are widely 
spaced, the spacing being dependent on the distance from the laser sensor and 
the increments of pointing angles used. Also because of measurement noise in 
the range readings there occur some false detections at few points. From 
Figure 3, it is clear that the outer edges will be detected easily, however, 
some of the inner edges may be missed for which the change of slope is small 
or gradual. Thus R.E.S. alone is not completely sufficient for extracting edges 
of the objects. 

Having described the surfaces in the scene by the least squares regression 
planes, the continuous and exact inner edges can be obtained by calculating 
their lines of intersections. It is to be noted that this approach can be 
applied only for the inner edges where both the intersecting planes are visible 
to the laser sensor. For outer edges the best estimate is a straight line 
passing through the edge points lying in a three-dimensional space. 

In Figure 4, the plane has been defined uniquely by its unit normal 
vector n^ = pT + qj + rk and its distance from the origin, shown as d. T, j, 
k are the unit vectors and p, q, r are the components of n^ along these directions. 


The equation 

of the plane is written as n^-r. = d + px + qy + rz 


I*! “ ^ 

+ px + qy + rz 

(12) 

which can be 

put in the form 



r ^ r ^ 

(13) 

= a + 

bx + cy 

(14) 

where z = ^ 
r 

K - ^ ® 

>. h--^,c-- — 




Fig. 3a Laser beam following a flat surface 
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Fig. 5 


Flow diagram of the 3-Dimensional Scene Analysis System 
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In Equation 14, parameter a * ^ is a measure of the distance of the plane P 
from the origin. The equations of the planes obtained from the least squares 

A ^ 

method of surface estimation are in the form of triples [a, b, c]. 

Thus two planes [a^ , b^ , c-|] and [agi b 2 » C 2 ] will intersect if 

b« T 


°1 *^1 

T-- = ~ indicates that the two planes in question are parallel to each other 
and hence do not intersect. All the inner edges can, therefore, be obtained by 
the above procedure after the equations of the two planes forming it have been 
obtained. 

2. Object Reconstruction and Recoqnition 

In a paper on two-dimensional scene analysis [14], Guzman used heuristic 
rules to decompose a line drawing into objects. These heuristic rules were, based 
on the types of vertices which occur in a line drawing. 

Here we present a heuristic scheme for object reconstruction and formation 
based on input data containing depth information. This scheme will reconstruct 
plane faceted objects from a workscene described as edges, faces, and vertices 
in cartesian coordinates. 

a. Object Reconstruction Usinq Heuristic Rules 

Planes will be collected into the same object if they satisfy these rules: 

1. Intersect at a convex edge. 

2. Contain col inear half edges which have adjacent coplanar faces. 

3. Intersect at a concave edge where one of the vertices of the edge 
is on the outside edge and the concave edge is not col inear with 
the outside edge. 


These heuristics are pictorially described in Figure 6, where 
indicates that the two faces connectf are grouped into the same object. Figure 
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Figure 6 
Heuristic Rules 

A. Faces linked because they intersect at a 
convex edge 

B. Faces linked due to coplanarity 

C. Faces lined because they intersect at a 

concave edge where one of the edges vertices 
fall on an outside boundary ' 




17 


6A shows faces which are grouped because of the convex edge between them. 

Figure 6B shows faces which are grouped because they contain col inear edges 
with adjacent coplanar faces. Figure 6C shows an example of grouping according 
to Rule 3. 

Rule 3 is included in the set of rules because a laser rangefinder, which 
measures only depth, will not see an edge between objects which are aligned. In 
Figure 7, Part C will be grouped into the main body when viewed from the right 
because it contains part of face 1, and face 1 has convex edges which link it 
to the rest of the main body. Face 1 contains vertices A, B, C, D, E, and F. 

If this object is viewed from the left, we want the same grouping to occur so 
that the description of the scene will be consistent. When viewed from the left 
Part C can be grouped into the main body by Rule 3. Part A exhibits the same 
characteristic. Part B is not grouped with the main body since the concave 
edge is colinear with the outside edge of the main body. 

Using the real world as an example, objects on a table are considered 
separate objects unless they are aligned with an edge of the table, 
b. Convexity of Edges 

One of the most important features of the object reconstruction scheme is 
determining whether or not an edge is convex or concave. 

To determine if an edge is convex or concave, we define the unit normal 

/\ /V 

vector and for the edge in both planes intersecting at the edge. Figure 8. 
These unit normal vectors lie in the planes forming the edge, are perpendicular 
to the edge, and point to the inside of the plane. Then we form the unit vector 
pointing from the rangefinder eye to the midpoint for the edge in question and 

A A 

call this vector R. If the dot product of R and the sum of and Lg 

R • (L^ + Lg) 
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Part B 


Part C 



Figure 7 
Rule No. 3 




Figu 

Determinin 

If a is less 
is convex. 





20 


Is greater than zero, a convex edge 1s Indicated. If this dot product Is less 
than zero, a concave edge Is Indicated. Examples of how convexity Is determined 
are shows In Figures 9, 10, and 11, where Is the projection of R Into the 
plane formed by and Lg. 

c. Collnearlty and Coplanarlty 

All half edges are searched to see If any apirs are col Inear. For two half 
edges to be col Inear, three of the possible six unit vectors formed by pairs of 
the four vertices must be equal (Figure 8). We choose to use the pairs (A,B), 
(A,C), and (A,D) and say If: 


2 - ABS 


AB» AC 

TTMTFTWf 


-ABS 


ABoAD 

ffABU ifAD'l I 


< K 


(16) 


where K Is a small number def!<;ndent on the amount of noise, then the half edges 
are coll near. 

Half faces which are coplanar will have equal unit normal vectors cssoclated 
with the col Inear half edges. 

If l!L^-L2ll<K2 (17) 


where K 2 Is another constant based on measurement noise, then the half faces 
are coplanar. 

If Rule 2 Is satisfied, then the half edges are merged and the adjacent 

A A 

coplanar half faces are merged. In Figure 12, * L 2 Indicating that the half 

faces are coplanar. Also AB Is col Inear to CE, and both B and C are virtual 
vertices, satisfying Rule 2. Hence, F2 will be merged with FI and AB will be 
combined with CO to form AD. 
d. Results 

Computer simulation was performed on an IBM 3033 computer In a program 


written In LISP. A sample workscene Is pictured In Figure 13, where iff 



A 


R f > 0 R • L 2 > 0 

R #(L^ + 4 ) > 0 
edge AB is convex 



Acute Convexity 



R»l^>0 R»L2>0 

A /K A A 

IR • L^l > |R • L 2 I 
R e(L^ + 4 ) > 0 
edge A 6 Is convex 



Figure 10 Obtuse Convoxity 
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signifies the position of the laser rangefinder eye, which is the origin of the 
cartesian coordinate system in which the input data is recorded. 

The input is in the form of lists of features. A face is defined as the 
list of edges which bound it; an edge is defined as a pair of vertices and con- 
tains the information on the unit vectors, and a vertex is defined by this list 
of its coordinates. For example, face FI is defined as (EF DF CD CE) in the 
sample workscene. Also, edge DF is defined as ((FI 0 -1 0) (F2 0 0 1)) where 
(1 -1 0) is the direction of the unit normal vector associated with FI at edge 
DF. Vertex A is defined as (-6 -20 12), its coordinates relative to the range- 
finder eye location. Which vertices are virtual is determined using the fact 
that virtual vertices are the endpoints of only one edge while other vertices 
are formed by the intersection of more than one edge. 

The output contains the description of the scene upon completion of the 
object reconstruction scheme. 

The output was: 

VERTUAL VERTICES: (L S T M R Q) 

NEWEDGES: (HN BO AP) 

CONCAVE: (CD DH) 

CONVEX: (DE KV BO) 

OBJECTl: (F1 F2 F3 F4) 

0BJECT2: (F5 F6) 

The scheme actually worked as follows: 

•Virtual vertices were found 

•Col inear half edges were fixed creating NEWEDGES HN BO AP and merging 
faces F7 and F8 into faces F3 and F4 respectively. 

•FI and F2 were linked across convex edge FD 

•F3 and F4 were linked across convex edge BO 

•F5 and F6 were linked across convex edge KV 

•FI and F3 were both linked with F3 based on Rule 3 at edges CD and DH 



No other links were possible so that the groups of faces were then considered 
objects. Note that F7 and F8 do not appear in any object list since they have 
been merged away, 
e. Conclusions 

A heuristic scheme for three-dimensional object reconstruction has been 
presented above which used heuristic rules based solely on three-dimensional 
geometric considerations. We found that given three-dimensional input describing 
a work scene, the objects contained in the scene can be reconstructed using these 
heuristic rules. In cases where enough information is available, parts of the 
scene which are masked by occlusions can also be rebuilt. The reconstruction 
is consistent and allows for error due to measurement noise. 

TASK B. POLYNOMIAL SPLINE APPRAQCH TO TWO-DIMENSIONAL SIGNAL PROCESSING 

A state space appraoch to 2-D data processing using spline function is 
reported in the following pages. The state space notation of splines not only 
keeps the notation simple, compact, and uniform, but also it provides the crucial 
point in the recursive data processing. The fundamental principles and limita- 
tions on recursive smoothing splines are explored, which are neglected in the 
previous research. One special case which results in forward on-line recursive 
algorithm is presented. 

1 . The State Space Approach to 2-D Vector Processing Using Splines 

The 2-D spline smoothing problem is defined as follows. A set of noise 
corrupted 2-D data Z = {z. .; i * 1,2,— M, j = 1,2, — ,N} is given, we wish to 

* J 

estimate the original function values and its derivatives from these discrete 


noisy measurements. 

The measurement process is described by 


i = 1, — ,M 
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where Hj) is the sample value of original function which is assumed of 
(mi ed partial) derivatives up to certain order. is random noise with zero 
mean and finite variance R. .. 

* J 

The approximating function adopted here is the 2-D spline function which 
is defined in wide sense, a piecewise function of continuous derivatives up to 
certain order and each piece is defined only in one grid, for example, the regular 
square grid as shown in Figure 14. 

The reason for using spline function is that the approximating error 
analysis has been well developed in the numerical mathematics [9, 10]. Another 
advantage of this approach is that under certain conditions the recursive compu- 
tation can be realized, which makes the real time operation possible. 

The ent.re spline function in domain R as dictated in Figure 14 would be 


Soo(C.n) 

Syo(C.n) 






Ko < K 

Cx £ i £ Cz » 
^i 1 ^ 1 ^i+1 • 
^M-1 - ^ ’ 


ho £ n £ hi 
ho £ n £ hi 
nj £ h £ hj^i 


^N-1 1 ^ 1 ^f! 


The objective function which we wish to minimize is 
M N 


J = 


if] jf] "^ij ^^ij ’ 

^j+1 


M-1 N-1 

+ I 2 O-P 4 J 


"^i.l 


i»0 j=0 


Mj 


5- -'h^ 

1 0 

where Pj^'s are weighting parameters and 

0 < p,j < I. 


C(s)dCdh 


09 ) 


( 20 ) 


C(s) is the smoothness measure function usually defined as function of 
the high order derivatives of S(^,n). 
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1 

The first term in Equation (3) represents the noise effect at the grid points 
and the second term is a measure of smoothness of the approximating function. 

The selection of C(s) and S(C,n) could be independent, but the optimum S(5,n) 
which minimizes a special //C(s)dydX is often unique. 

Notice that S(C.n) would act as an interpolation function when p.. approaches 

' J 

to zero, and on the other extreme, when p.. approaches to one S(5»n) would become 

' w 

a plane. Thus, the value of weighting parameter p. • controls the trade-off be- 

I J 

tween the measurement error and the smoothness of the approximation function, 
a . The State Space Model of Bicubic Spline 

The most often used splines are polynominal splines because of its simplicity 

in computation and the closed form solution. In general, a 2-D spline defined 

in one grid is described by a partial differential equation. However, only few 

of the solutions satisfying some restrictions could be represented by linear 

discrete state space model defined on the grid points. The bicubic polynominal 

spline is the one among them of particular importance. 

The bicubic polynominal spline is the polynominal which minimizes the 

r 3^s T 2 

smoothness measure C(s) = ~~ — s- . For convenience, new notations of 

[kW} 

variables are introduced 



^i 

- ^ i ^i+1 

X A n - f 

"j 

< n < 

thus. 



0 £ u £ AC , 

AC 

A^i.r^i 

0 < X < An , 

An 

A 


From the principle of calculus of variation, the optimum solution of 



30 


Is the one satisfies the Euler equation, 


3®S..(y.X) 

L3 = 0 


with proper boundary conditions [15]. 

If only polynominal solutions are of concern, then the solution satisfies 
Equation (22) is equivalent to the one which satisfies the following equation. 


( 22 ) 


3'*S..(y,X) 

^ = [yX X y 1 ] 


3y^3X2 




(23) 


where ur^ ulj are under-determined parameters. 

The general solution of Equation (23) is 
S^.j(y,X) * (}).|(X) + yc()2(X) + 'i'-](y) + X((;2(y) 




y^X^ y^X^ y^X* 1 

4 J 


(24) 


'ij 


where (j)^, (f> 2 , il>y andij^ are aribitrary polyncmiinal functions. If we define the 

"state" as the column vector [S, . ||- , exact linear 

discrete model of Equation (24) defined on the end points of one grid can be 
obtained. 




(25) 


where 


s CS(y.A) (y.X) If (y.X) |^Jy,A) 


Mn,j+m ^ 3y 3X MX'"'*''' ■' 'XwnAn 


(25a) 


t « 0, 1 ; m * 0,1 
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poi = 
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0 0 
0 0 


0 0 
0 0 
1 A, 
0 1 
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0 
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1 

0 


0 

s 

0 

0 

1 , 


(25b) 


and 


:0 0 - 


r = 


Aj-A 
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12 

12 
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V? 

A^A^ 
C n 


12 

6 

4 

2 

A^A^ 

AtA^ 
%. n 

A^A 
C n 


12 

4 

6 

2 

A|A^ 
L 4 

A-A^ 
C n 

2 

C n 
2 

A.A 
? n 


(25c) 


(25d) 


The linear discrete state model for a general 2-D spline Is often not 
available. Because of the restriction of the positions fo grid points, it only 
allows the function S(y,A) propogates the variation along the lines parallel to 
the C or n axis. Consequently, the arbitrary functions could only be functions 
of only one variable, y or X, so that they can be cancelled. Linearization or 
approximation may be used to derive the linear discrete model for those do not 
meet above requirements. 
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Since the finite state model is one of the essential condition for 
recursive algorithm, the bicubic spline is the one, if not the only one, of the 
members which can be implemented recursively. 

b. Vector Processing Formulation 

1-D estimation theory is well developed, nevertheless, the corresponding 
theorems in 2-D is not always possible. This complication occurs partially due 
to the choice of the support defined for the 2-D filters [16-18]. An estimator 
with support defined by Equation (25) usually does not have the optimum 
solution either in the sense of least square or quadratic performance index. 

To avoid this fundamental restriction, a vector processing model is 
constructed as follows. Define the glcbal state quantity to be the vector 
along the column of region R, 






(26) 


j=0,l,2....N 





The state equations along the column strip are 

" rUH,j 




(27) 


or 




+ F^^X. 
" -J 


+ F®®X. 



(27a) 
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which could be rewritten as 

ij+l ■ f ij ♦ r Uj ( 28 ) 

where 

F « [I - [F^® + F°®] , 

r » Cl - r® , 
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To simplify the calculation and deduce the recursive solution, only the 


quadratic smoothness measure is considered, that Is from Equation (6) assume 


0 JO 


rA fAj. d 

^ C(s)dydA * ^ 

Jn Jn Jf 


3i?wJ 


dudA 


’i] 8,j 


where A 


rA, 


A 


[pA A p 1 [pA A p 1 ]dpdA 


(29) 

(29a) 

(29b) 


for bicubic splines. However, the result of above problem can be easily extended 


to the general quadratic form. 


C(s)dpdA.[U^. X.. 


"ij 

^1+1 J 
X1,j+1 


(30) 
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r 


Let zj A [0 2 ^j Zgj • 
then the objective function (Equation 20) can be expressed In the global state 
form, 


■ <*o ■ ic> ^o’ I £J' 

j*l 


N 


iT „-l 



where 2^ Is the initial guess. 


H = 


H 0 O 

OH. 0 

* 

. O o’ ■ H 




° o 

''tj “ 


0 O-»,j)0,-] . 

0 O 


O 

0 




(31) 


Ola) 


(31b) 


(31c) 


The first term In Equation (14) Is added to reduce the effect of initial guess, 
c. Smoothing Algorithms 

We now turn to the development of the estimation algorithms. The closed 
form recursive formula for the problem described in the previous sections does 
not always exist. There are some significant special cases that provide simple 
recursive computations as discussed below. 


35 


For the sake of forward real time recursive computation, the unknow 
coefficients Uj can only be the function of previous states J^^'s and parameters 
y^'s, i<j. Therefore, the minimization procedure of Equation (31) can be 


expressed as 
Min 


jJ . Min 
X , 

-0 -0 — J -1 


Jo* 


IHn 


J, * 



(32) 


where Z'j A {z^^: i 1 1 £ j } 


and 


\ k (^t - H (h 


a Xj) ♦ il[., aj., 

I “ 1 ,2, j 


A ^ AT 1 ^ 

‘(]L‘ Lr P ' (X - x^) 
0 -0 -0 -0 0 


j 


(32a) 

(32b) 

(32c) 


1. 


t=i 


The notation |Z'j emphasizes that the measurements are up to stage. 


Let us denote the optimum estimate of which minimize is 
i*» li| j is the optimum estiamte of using the data z^, Z 2 tZj. The 

* A A 

relation between and j comes from Equation (25) 


iitiij 




^^113 


where the index i/j indicates the quantities are the optimum solution using Z^ 


(33) 

J 


In order to obtain Xjjj, we must compute all U^jj's and X^^ via Equation 
(32). The whole procedure of computation is dictated in Figure 15. First, we 
compute Uj|j's from stage backward to 0 ^^ stage, and then the smoothed 
estimate S^jj's are calculated forward by using Equation (33). 

At the state, Uj_ijj is obtained by minimizing Jj with respect to 
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'j 


= (ij -HXj)'£:’ fij -HXj) Uj., 


T „-l 


=(Zj-HFXj., -HrUj_,)'R-‘ (ij-HFXj.,- 

*M]., 2]1, 




( 34 ) 


since V 


, ^ ^j -1 ^j-T taking the gradient of Equation (17), we 

obtain the expression, 

-i'il r:'(zj - £ f X.., - H r U..,) . ajl, Uj., = o 


or 


Mj., = T-’r R-’ (Zj -HFXj,,) 


Where!.- £ 

sj vl * yJ 


(35) 

(36) 
(36a) 


Because the proceeding solution is solved by using data Z'? , U._,, and X. , should 

• J * • J ” • 

be denoted as and respectively. 

To compute the U._p Equation (36) is substituted into Equation (34), and 
J. is simplified as function of , only. 

J w ' 


Jj • (£j -HFXj.,)Tw-> (z. - HFXj.,) 
where + H ! £j_., ^ 

At stage, the problem to be solved is: 


(37) 

(37a) 


mm 

!ij-2 




• mn [Zj_, - H Xj ,)’^ Rjl, (Zj_, -HXj_,) 

- 3-2 

a]!2Mj.2MZj -HFXj_,)^W-' (z.-HFXj.,)] 


The resultant Uj_ 2 |j is 


(38) 


2 j.2jj = Ijll 'l' / sjll (ij-l - H F Xj.j) * rT f «:'(Zj - H F FXj.^)) 

(39) 
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where L.i = {iVsi!? M £ (39a) 

Observing Equation (36) and (39), it seems that there is no simple closed recursive 
formula to compute 

However, if the state Equation (33) satisfies certain conditions, the closed 

/N 

recursive formulas for U. j^'s or^^. j^.'s are possible- Substituting Equation (22) 

to Equation (21), becomes function of Xj _2 only. Repeat the above procedure 
th i 

backward down to 0 stage, we will have which is a function only. Then 
>^jj can be obtained by minimizing J“ with respect to After that, X-]|j, 

X 2 jj, ... are then computed forward via Equation (33). This nonrecursive computing 
procedure can be applied to the general state model on which no restrictions are 
posed. 

There are two special cases which save the computation and/or memory storage 
significantly. The first case is that the calculation of jj could be recursive, 
if H is invertable. However, in many practical applications, this may not be 
true. 

The second case is that we can compute X^. jj directly from without 

computing all the U^-jj's, i= 0, 1,..., j-2, if Pis invertable. The second case 
is much more important not only because it is true for the splines defined in 
Section a , but also it r sjlts in a true forward recursive on-line algorithm. 

Assume r is invertable, U^. * r’^ (£^^^ - FX^. ). Because all IJ^. 's are inde- 
pendent variables, the jcps can be considered as independent variables. There- 
fore, J^. can be rewritten of the form 

Ji = (li - I X.)^R'^ (z. - H _X.) 

+ (ii - ^ I ^ (X. - X._^) 


(40) 
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Thus, 


min 


Oq = mi n 


^0, 


= min 


ij-Aj., 


J. + min 

J 


h-uk-2 


j-1 .+ + min J. 


C41) 


'0 


Noted that every term in Equation (41) is a function of two variables only, and 
every indpendent variables 's are involved in two terms only. In comparison 
with Equation (32), Equation (41) has a significant decoupling structure which 
enables the direct forward recursive algorithm. 

To illustrate this feature, let us see the following example. Assume the 


optimum filter estimate Xjij is already available, which is the solution of 


min 
^0 * 


mm 


J . + mi n 

J 


Jj _1 + + min Jq 


(42) 


*j-rV2 


Now, we wish to estimate , that is, we want to find the X..j.i such that 


min 

^0 * x^ , . . . . Xj ,x j^1 


= min 




J. , + min 

J+1 


Compare Equations (42) and (43), it is evident that 


J. 


+ mi n J, 


(43) 


min 

^O’^l ^j+1 




= min 
^j+i 


J + mi n 


iJ 

'^0 


(44) 


^0 *^^ ****** ^ J 

Consequently, the present estimate would equal to the previous estimate 

/\ 

Xj|j plus a correction term which depends on only. 

After some derivation, filtered estimate is given fay 

hTr-:, (ij,, -HFX.|j^ 


(45) 



«J.1 = £ 2j J * £ £j 

Similarly, the 1-lag smoothing equation can be obtained, 


-jIj+i “-jlj '^Ai= = =j+l ^-j+1 "=-j+l|j+l^ 

which can be extended to m-lag smoothing equation, 

h = X. . + p. h'*’ R:^ (z. - H X. 

-*j-m|j — j-m j-m =j-m = *= *=j-m+l — j-m+1 =— j-m+1 j 


(45a) 

(45b) 

(46) 

(47) 


d. Conclusions 

The general properties of 2-D smoothing splines problem has been discussed 
in this paper. The principles and restrictions on forming a recursive algorithm 
of 2-D vector processor using splines are f-x'l-tred. It is found that the recursive 
on-line structure for general splines is not always possible. One special case 
which result in forward on-line recursive algorithm was developed. 

The adoptation of state space model clarifies the crucial point of recursive 
spline smoothing algorithm, which were neglected before. However, the adjustment 
of weighting parameter p. . is not included in this report. The development of 
adaptive algorithms and the extension to nonregular grid splines will be presented 
in the future. 


2. Basis Functions For Piecewise Hermite Polynomials (p.H.p) 

A systematic method for deriving the bases of piecewise Hermite polynomials 
is reported here which can be extended for the general cases of 2-D and higher 
dimensional p.H.p. 's. Among the polynomial splines, piecewise Hermite poly- 
nomials are of particular importance because they have the following attractive 
properties: (i) The p.H.p. requires data in only a finite local area around the 
point of approximation, hence it can be implemented recursively (finite local 
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support), (ii) The continuity of p.H.p. can be prespecified which results in 

a smooth function for curve fitting, (iii) p.H.p. is of finite error bounds and 

the error behavior can be i lalyzed accurately. We denote the collection 

n n+m 

cf all piecewise 2-D functions that have continuity up to nth order partial 
derivatives on the first variable (c)i mth order partial derivatives on the 
second variable (n) , and the n+mth mixed partial derivatives on these two variables. 
Because of the requirements of continuities on the boundaries of every piece (grid), 
the order for the polynomial which meets these boundary conditions are 2(n+l) 
for the first variable and 2(m+l) for the second variable, which is the so-called 
"piecewise Hermite polynomial (p.H.p.)." 
a. Cardinal Method 

Consequently, the p.H.p. of order 2(n+l) of polynomials of the first variable 
and order 2(m+l) of the second variable would form a basis of the linear space 
n*'n+m‘ total number of the p.H.p. bases is 4(n+l )(nt+-l ) in each grid. 

If a 2-D partition n(C,n) called observation region is given as shown in 
i-igure 16, with (N-1) intervals in the variable and (M-1) intervals in the 
variable n, the p.H.p. of each grid defined by 

^ 1 ^i+1 » ’Ij i i ... (N-1), j=l,2, ... (M-1) 

can be constructed by the following procedure. 

(i) Because of the independence of the first and the second variables, the 
p.H.p. basis can be formed by the tensor produce [6] of the 2(n+l) p.H.p. basis 
factors (polynomials of degree, 2n+l ) in the first variable C, (5). and 
the 2(m+l) p.H.p. basis factors in the second variable n, (n). Hence, 

the general p.H.p. basis in one grid can be expressed by the combination of 
(5) and (n). 
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1 

(ii) According to the continuity requirements along the boundaries (4 end 
points of the grid), a set of 4(n+l) (m+1) equations are posed on these p.H.p. 
basis factors. 

(iii) Each of the p.H.p. basis factors is obtained by using the boundary 
and normalization conditions. 

b. Bases for Twice Differentiable Continuous Functions in Both Directions 

To illustrate the above procedure, an example for consturcting the p.H.p. 

2 

bases of is given in the sequel. 

If a function f(C,n) is constructed by using the p.H.p. bases from the given 
data vector 

.33 3^f 3^f 3^f 3^f 3^f 3^f _ „T 

3C * an ’ an ’ ’^2 ’ ^,2 ’ 3^2 3 , ’ 35 3,2 ’ 3^2 an 2 " " 

at every grid end point, the resultant function f(C,n) will satisfy the continuity 

2 2 
requirement of and thus belong to 2 ^^. 

(i) The whole function can be expressed as the summation of piecewise 

functions in every grid, 

N-1 M-1 

f(C,n) = E I ^ ^ (^^8) 

» »J 

i -1 j=l 



( 49 ) 
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where 

♦ok(5)‘'t'2ii(5) • 't'2k*5>'*'U*'’) ■ ‘•'lk*^>'*2ltf'') ’ ♦2kf^*''*’2lt''*^ 

4)0k(O. 4^k(C), p.H.p. basis factors of 

variables C and n. respectively. Each of these is a polynomial of degree 5. The 
complete p.H.p. bases are shown in Equation (50). 

The next step is to determine every p.H.p. basis factor mentioned above. 

Each p.H.p. basis factor is a fifth degree polynomial, hence, there are six 
undetermined coefficients for each factor and a total of 72 parameters to be 
calculated for these 12 basis factors. 

(ii) From the given boundary conditions, these coefficients are determined. 
This results in 12 sets of simultaneous equations corresponding to 12 factors 

where D? = d*^/d^*^ and = d*^/dn^ , <S. u "'S Kronecker 6 function 

c, ri a ,D 

p = 0,1 ,2 ; q = 0,1 ,2 ; and k = 0,1 ; d = 0,1 . 

(iii) By solving the simultaneous Equations (51) and (52), every undetermined 
coefficient can be obtained. Because only one term is nonzero for each basis 
factor (for a fixed q and k in Equations (51) and (52)), it is easy to obtain all 
the six coefficients of this particular fifth degree polynomial. 

It can be written without difficulty that 

♦oo(c) = 
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The unknown parameters aQ, bg, and Cq can be calculated by the equalities 
4>oo(?i) * ° simplicity, let A? A - C,- = 1 

and X * $- then the foregoing equation becomes 

(j>Oo(^) = (X-1)^ (-6X^-3X-1) = 1-10X^+15X^-6X® 

By the same token, all the other basis factors can be derived similarly. 


= X - 6X^ + 8X^ - 3X^ 


. m - 1^2 3,3 3,4 1,5 

<(>Q^(X) = lOX^ -15X^ + 6x^ 
<j)^^(X) * - 4X^ + 7X^ - 3X^ 

4»2^(X) ^X^ - X^ + ^X^ 


(j;Oo(v) = 1 - lOv^ + 15v^ - 6v® 

ij^io(v) = V - 6v^ + 8v^ - 3v^ 

(pQiC'^) = lOv^ - 15v^ + 6v^ 

<1^1 l(v) = - 4v^ + 7v^ - 3v^ 

'i'2i(v) = 




4 1 5 


where An = n^_j_^ - n^ = 1 and v = n- n.,*. 

After obtaining all the 12 basis factors, the p.H.p. bases are readily 
obtained by combining these factors as indicated in Equation (50). 


c. Conclusions 

A systematic method for deriving the bases of p.H.p. 's has been described. 

This approach can be used either for high order polynomials or high dimensional 

functions. As can be observed from the example, the derivation procedure and the 

calculations are simple. Consequently, this procedure is suitable for practical 

2 

applications. For example, the piecewise Hermite polynomials of can be used 
to estimate the second derivations of 2-D image. Because of the smoothness pro- 
perty of piecewise Hermite polynomials, they can also be used to filter out the 
noise of a stochastic signal. 
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